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O 1 ABSTRACT 

CO 

We have studied the impact of cosmic-ray acceleration in SNR on the spectra 
of cosmic-ray nuclei in the Galaxy using a series expansion of the propagation 
equation, which allows us to use analytical solutions for part of the problem and 
an efficient numerical treatment of the remaining equations and thus accurately 
describes the cosmic-ray propagation on small scales around their sources in 
three spatial dimensions and time. We found strong variations of the cosmic-ray 
nuclei flux by typically 20% with occasional spikes of much higher amplitude, but 
only minor changes in the spectral distribution. The locally measured spectra 
of primary cosmic rays fit well into the obtained range of possible spectra. We 
further showed that the spectra of the secondary element Boron show almost no 
variations, so that the above findings also imply significant fluctuations of the 
Boron-to- Carbon ratio. Therefore the commonly used method of determining CR 
propagation parameters by fitting secondary-to-primary ratios appears flawed on 
account of the variations that these ratios would show throughout the Galaxy. 



Subject headings: cosmic ray propagation, cosmic ray origin 
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1. Introduction 

The origin of cosmic-rays (CR) is one of the major problems in modern astrophysics. 
In that quest two possible strategies may be followed: 

• one may search for high-energy emission from the source regions themselves, for the 
CR density must be very high therein. Particle acceleration at supernova remnant 
(SNR) shock waves is regarded as the most probable mechanism for providing Galactic 
cosmic rays at energies below 10 15 eV on account of the power requirements (Blandford 
& Eichler 1987). The recent detections of non-thermal X-ray synchrotron radiation 
from the four supernova remnants SN1006 (Koyama et al. 1995), RX J1713. 7-3946 
(Koyama et al. 1997; Slane et al. 1999), Cas A (Allen et al. 1997; Gotthelf et al. 2001), 
and RCW86 (Bamba, Koyama, and Tomida 2000; Borokowski et al. 2001; Rho et al. 
2002) support the hypothesis that at least CR electrons are accelerated predominantly 
in SNR. The evidence for CR nucleon acceleration in SNR is at best controversial 
(Berezhko, Ksenofontov, and Volk 2002; Reimer & Pohl 2002). 

• one may also construct models for the diffusive propagation of CRs in the Galaxy and 
compare their predictions with observed CR spectra at earth (e.g. Strong & Moskalenko 
1998; Maurin et al. 2001; Taillet & Maurin 2003). The relation between primary CRs, 
which have been accelerated somewhere in the Galaxy, and secondary CRs, which 
are produced in interactions of CRs with ambient gas, then allows to determine the 
propagation parameters such as the diffusion coefficient, the size of the diffusion region, 
and others. Nearly all published studies of that kind assume steady-state conditions. 

Already in 1969, Lingenfelter (1969) found the CR energy density at earth to vary due 
to nearby SN, given the CR are accelerated at these sites. A few years ago Pohl & Esposito 
(1998) showed for CR electrons that in case of SN origin the CR electron spectrum would 
vary in space and time throughout the Galaxy. The purpose of this paper is to present a 
time-dependent propagation model for CR nucleons. In particular we aim to address the 
following questions: 

• would a SNR origin of CR nucleons also lead to significant fluctuations of the CR 
density in the Galaxy, which then would modify secondary-to-primary ratios from 
their steady-state values? If that was the case, we would have to re-think the way we 
infer CR propagation parameters from their locally observed spectra. 

• are there any signatures in the CR distribution in the Galaxy, that might permit to 
infer a SNR origin of CR nucleons on the grounds of locally observed CR spectra and 
the diffuse Galactic gamma-ray emission? 
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For that purpose we have developed a method to numerically solve the propagation 
equation based on a series expansion, which allows us to use analytical solutions for part of 
the problem and an efficient numerical treatment of the remaining equations. As described in 
detail in section 2 this technique allows us to accurately follow the CR propagation on small 
scales around their sources in three spatial dimensions and time. In section 3 we discuss 
the results of these calculations with respect to the CR density distribution and in section 4 
with respect to CR spectra, respectively. 



2. The Model 

2.1. The basic equations 

Our investigations are performed in the framework of a diffusion model of CR propaga- 
tion, so we use a continuity equation for the differential CR number density, N: 

dN 

— - S = V (WiV) - QvaN. (1) 

with particle speed v and total spallation cross section a. We consider the diffusion zone to 
be a disk of radius R = 20kpc (cf. Webber et al. 1992) and height 2H (see Fig. 1), where H 
is determined below. The density of the interstellar gas is assumed to decrease with distance 
from the galactic plane, z, with a profile 

The parameters no = 1.24 cm -3 and /i 9 = 30kpc~ 1 corresponding to a column density of 
~ 6.2 • 10 20 cm~ 2 are chosen to be consistent with the calculations by, e.g., Webber et al. 
(1992) or Berezinskh (1990). The chemical composition of the ISM, which consists mainly 
of hydrogen and some heavier elements, is taken into account by multiplying the gas density 
Eq. (2) by a factor of 1.3, as it is done by Mannheim & Schlickeiser (1994) for the calculation 
of CR-induced pion production in the ISM. 

The diffusion coefficient, k, depends on the particle rigidity ( = p/q with particle mo- 
mentum p and charge q. The form 

ko (£) for C < Co 



with Co = 4GV/c is chosen to reproduce the observed Boron-to-Carbon ratio. 
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Fig. 1. — Sketch of the geometry used in our calculations. We assume that the Galactic 
disk with radius R is filled with interstellar gas and the CR sources. The density of the 
interstellar gas decreases quasi-exponentially in the halo. The height of the entire diffusion 
zone is 2H. 
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As we consider discrete, point-like sources, the source term S in Eq. (1), is in fact the 
sum of the contributions of many source, each of which has the same temporal and spectral 
form. For the time dependence we assume a linear growth with an exponential cut-off, with 
the source spectra being power laws with index s in particle rigidity (. 

Qi (C, t) = qj ■ (t - tj) exp (-^yt) © (t - tj) ■ (£) . (4) 

We assume the SN explosions to be stochastic events. Their effect on the CR spectra can 
thus be studied by the method of Monte-Carlo simulations, so many possible CR spectra 
are calculated using randomly chosen sets of CR sources. We used the rejection method 
(Press et al. 1993) to compute the quantities rj,{pj,Zj for the location, tj for the ignition 
time, and (jj for the source strength. The position of the CR sources in azimuth, ipj, is 
uniformly distributed, whereas for the radial distribution we use the form suggested by Case 
& Bhattacharya (1996), and for the vertical distribution we use the same profile as for the 
density of interstellar gas. Then 

P s (r,z) = 1 ( —Y- exp ( ar °- br ) . (5) 
cosh(z h g ) \ar s J \ r s J 

with parameters a = 1.69, b = 3.33, r s = 8.5 as taken from the paper of Case & Bhattacharya 
(1996). 

The ignition time, tj, and the source strength, qj, are uniformly distributed on the inter- 
vals [tstart; ^cnd] and [0; 1], respectively. We have also performed simulations using a detailed 
model of the spatial and temporal evolution of the nearby star-forming region Gould's belt 
(Perrot & Grenier 2003). 

This leaves k , the halo height H and the source spectral index s as free parameters, 
which we determined by fitting the Boron-to- Carbon ratio and the survival fraction of 10 Bc 
calculated for the steady-state case to measured data. Solar modulation is taken into account 
using the force field approximation (Gleeson & Axford 1968) with modulation parameter 
$ = 500 MV. In the steady-state case the best-fit parameters are k = 0.26 kpc 2 Myr~ 1 , H 
= 2kpc, and s = 2.1, as shown in Fig. 2. 



2.2. The method of computation 

To determine the influence of the discrete nature of SNR as CR sources on the CR 
distribution and spectra, one has to solve the CR propagation equation Eq. (1) with high 
resolution, both in space and time. With today's computers, however, it is not possible to 
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Fig. 2. — Best fits of Boron-to-Carbon data (upper panel) and 10 Be/ 9 Be data (lower panel), 
compared to data by Engelmann et al. (1990) (A), Dwyer & Meyer (1987) (O), Krombel & 
Wiedenbeck (1988) (□), and Orth et al. (1978) (x) for the Boron-to-Carbon ratio and from 
Connell (1998), Garcia-Munoz et al. (1981), Hams et al. (2001), Lukasiak et al. (1994), de 
Nolfo et al. (2001), and Wiedenbeck & Greiner (1980) for the 10 Be to 9 Be ratio, respectively. 
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calculate the CR density with a spatial resolution sufficient to investigate point-like sources in 
an adequate way. One should note that the grid size in finite-difference algorithms, e.g. those 
in the widely used code Galprop (Strong & Moskalenko 2001), in principle has to be much 
finer than the size of the presumed sources, for one has to deal with large gradients in the CR 
flux on small scales. Thus, the accuracy of representation of the CR density on a fixed grid is 
highly questionable. One way out of this quandary may be the implementation of adaptive 
mesh-refinement techniques in the finite-difference algorithms. The way we followed is to 
use an analytical ansatz that breaks down the problem of solving the propagation equation 
Eq. (1) into many small tasks, that are easily solved on pc-style hardware. The series ansatz 
to solve Eq. (1) will be in detail described in the next section. We also developed a method 
to spread the task of solving the propagation equation onto many computers, so the spatial 
resolution obtained is only limited by the number of computers available. 



2.3. The Ansatz 

Considering the cylindrical geometry of our Galaxy (with radius R and height 2H, see 
Fig. 1), we assume the gas distribution Q and diffusion coefficient k to be independent of r 
and ip. Then the CR transport equation Eq. (1) can be written in cylindrical coordinates: 

ON n ,,.[ldN d 2 N 1 d 2 N d 2 N) n . , Ar 

We start by expanding the desired solution N in a Fourier series in ip and a Fourier-Bessel 
series in r: 

N = -VV (A nm • cos (nip) + B nm • sin (nip)) Jn """f 2 (7) 

W nm {Jn( a nmR)) 

with a nm being the mth solution of j n (a nm R) = (in ascending order). We also expand the 
individual source terms for the point-like sources 

S = 5>(p, (8) 

i 

where qi(p,t) is the time and momentum dependence of the particular source and r-i its 
position, into the same series in r , ip 

S = ^iMSiz-Zi) (9) 

i 

x > > (cos (nip) cos (nipi) + sin (nip) sin (nip { )) — — — - 2 . 

(j'n (anmR)) 



n m 
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For extended and continuous sources, e.g. for secondary CRs, one has to use the Fourier- 
Bessel representation of the source distribution of the particle in question. 

Inserting Eq. (7) and Eq. (9) into Eq. (6) and using the orthonormality of sine and 
cosine and the analogous property of the Bessel functions (Watson 1944) 



R 



Jn {p^nmT) Jn \Ol n^T ) 



[j'n (&nmR)) 

one obtains equations for the expansion coefficients A 



rdr = 5 (m, /x) 



dt 



— ^ y 



Kp) { -^ 2 nm A nm 



d 2 K 



dz 2 



— Q(z)v aA r 



with 



cos [n ifi 



{j'n (oinmO))' 



and similar equations for B„ 
dB ri 



dt 



— 



HP) i -almBnm 



+ 



d 2 B r , 



dz 2 



— Q(z)vaB r , 



with 



S 



B 



E 



qi(p,t) sin (nipi) 



Jn i^^nm^i) 
(j'n (ttnmO)) 1 



(10) 



(12) 



(13) 



(14) 



These equations are not analytically solvable for arbitrary Q(z), therefore one has to use 
numerical methods. The advantage of this ansatz is, that it is possible to solve these one- 
dimensional PDEs simultaneously on many computers. Also, the resolution at a given point 
obtained in r,ip merely depends on the number of coefficients in the series ansatz, Eq. (7), 
that is actually computed. 

As the aw's increase monotonically with n and also with m, one sees from Eq. (11), 
that for large m,n we have k{p)a 2 ^> Q(z)vaB nm and therefore, the latter term may be 
neglected. In this case an analytical solution is known. 



2.4. CR distribution due to a single source 



Before we start solving Eq. (11) and Eq. (13), we have to determine the number of 
coefficients to be taken into account in the ansatz Eq. (7). Evidently, the spatial resolution 
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in r, <f direction of the solution obtained by the ansatz Eq. (7) depends on the number of 
coefficients used and also on the distance r from the origin. As we are mainly interested in the 
CR density in the vicinity of the Sun, we determine the number of coefficients by comparing 
the solutions of the propagation equation Eq. (6) obtained by a truncated series according 
to ansatz Eq. (7) with the solution of the propagation equation for a spherical source with 
a radius of p s = 50 pc, located at the position of the Sun. To ease the calculations we 
neglect the geometry of the Galaxy, assuming the loss processes not to depend on the spatial 
coordinates. This is permitted if we consider only the vicinity of the source for a limited 
time after the SN explosion. So, we study the source in an uniformly distributed ISM. 

First, we derive a solution for a spherical source. In this case, placing the source at 
the origin of our coordinate system provides us with spherical symmetry, i.e. the solution 
only depends on the radial coordinate, p. We further restrict ourselves to one fixed particle 
momentum p. Then Eq. (1) can be written as 

ON n , fd 2 N 2 8N\ , Ar 

in- s = k W + - P W- m (15) 

with b — flva. The source S has the temporal form given in Eq. (4), so we have 

S = M*-^)exp(-— W*-^)e>(P*-p)- (16) 



T 

Without loss of generality, we set qi = 1 and tj = 0. 
For Eq. (15) one can find the Green's function 

G = ) exp(-6 (t- t )) 

8ny/nVk ppoVt - to 

For a verification that Eq. (17) is indeed a Green's function for Eq. (15) and a short discussion 
on how to obtain it, we refer to Appendix A. 

Thus, we obtain the solution of Eq. (15) with sources Eq. (16) by the convolution 
N(p,t) = f [ G(p,po,t,t )-S(po,t ) dp o dt (18) 

Jp =0 Jt =0 

= r f 7 ex P (~ b (* - to)) ex P (~ ; f! + ^ x 

Jp =QJt =Q%iry/Try/kppoy/t-to \ 4k(t-t ) 



XSinh (^F^))' t0eXP ("r 

xe(t )e(p s -p)p 2 dp dt . (19) 
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The po integration can be performed analytically, which leads to 
N(p,t) = f (-2 y/k{t-t )exp( ' ~ ps2 ~ p2 

Ps- P 




Ak (t-t ) 




2 y/k{t-t ) \2 y/k(t-t ) 




xexp (b (t - t))t exp ( — — )dt . (20) 



T 



Unfortunately, this integral is not solvable analytically. It was computed numerically using a 
Riemann sum. The result is plotted in Fig. 3, where we also plotted the solution of Eq. (1), 
using addends with m < 210 and n < 311 in ansatz Eq. (7), which turned out to be the best 
tradeoff of numerical complexity versus spatial accuracy. This solution using ansatz Eq. (7) 
describes well, both in r— and (^-direction, the spatial evolution of the CR density around a 
spherical source. 



3. Signatures of Discrete Sources in the Galactic Cosmic Ray Distribution 

In this section we represent first results of our investigation as to what extent the SN 
origin of CR affects the density distribution of Galactic CR. We study the temporal evolution 
and the fluctuations in the CR density, using the methods developed in the previous section. 



3.1. Randomly Distributed Sources 

We now consider the time-dependent propagation equation Eq. (6) that we solve using 
the ansatz Eq. (7) for several source distributions. As a first step we performed a calculation 
for CR point sources that are randomly distributed in the Galactic plane. For that purpose 
we chose 16 O, for it is regarded as the most abundant primary CR element beyond helium. 
For this calculation and those described in the following sections, the CR density N was 
computed using ansatz Eq. (7). The corresponding coefficients A nm , B nm were obtained 
numerically using a semi-implicit scheme based on the Du Fort-Frankel (Du Fort & Frankel 
1953) and leapfrog schemes. We started the calculation of the temporal evolution of the 
CR density ,N, from the appropriate steady state solutions, which were calculated using the 
package TOMS638 (Houstis et al. 1985a,b). 

For the source function we use n q randomly distributed point sources with a radial 
probability distribution given by Eq. (5). Considering only one CR primary element at a 
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Fig. 3. — Comparison of CR density due to a spherical source with radius 50 pc given by 
Eq. (20) (solid line) and that computed using Eq. (7) with 312 coefficients used for the series 
in n and 210 coefficients used in the series in m; cuts in r (dotted line) and ip (dashed line) 
direction, for different times from source "ignition" . 
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time, we have for the source term S in Eq. 



(6): 



S(r,ip,z,t) 



^2qi(t,p)5(f-fi) 



(21) 



i=i 



with the source strengths defined in Eq. (4), and r*j the locations of the ith source. The 
number of sources, n q = 150000, was chosen to reproduce the local supernova surface density 
of 25 SN per Myr and kpc 2 of the Galactic disk (Grenier 2000). We calculated CR densities 
for several random SN distributions in the Galactic disk, each for a time interval of 10 Myr. 

In these calculations, the coefficients with 1 < m < 210 and < n < 311 of the series 
ansatz Eq. (9) were computed. This yields a resolution in r, (p of ~ 170 pc at the position 
of the Sun. The grid spacing in z is 20 pc, the time step is 1 kyr. 

Computing the series Eq. (9) at the position of the Sun gives the time variation of the 
CR density at this position, as is shown in Fig. 4, where the CR density for an energy of 
10 GeV per nucleon (upper panel) and for an energy of 5 TeV per nucleon (lower panel) 
is plotted versus time. These figures illustrate that the density of CR with an energy of 
5 TeV shows more rapid fluctuations than at an energy of 10 GeV on account of the energy 
dependence of the diffusion coefficient. 

The variations in the CR density at a given location, shown in Fig. 4, have a typical 
amplitude of about 20 per cent with occasional spikes reaching 100 per cent. The latter 
mainly occur due to nearby SN explosions, as illustrated in Fig. 5, where the temporal 
evolution of the CR density at an energy of 10 GeV per nucleon is shown for a 400 pc x 400 pc 
section of the Galactic plane. Variations of the same order of magnitude of the local CR 
energy density due to nearby SN have been found by Lingenfelter (1969) when he investigated 
the contribution of these SN to the CR energy density near our Sun, using a simple diffusion 
type model only considering losses due to escape. The perturbation in the CR density due 
to the source stays almost localized and does not spread out. 



The local distribution of stars and the ISM is dominated by the so-called Gould's Belt 
(Poppel 1997). Of particular interest for our investigations is the enhanced SN rate within 
the belt (Grenier 2000). The kinematics of the belt can be modeled by an expanding ring 
of gas, assuming an initial explosive event (Perrot & Grenier 2003). Its age is estimated to 
be about 30 Myr which implies that with regard to our calculations, which cover the last 
10 Myr, not only its presence, but also the still evolving geometry of the belt has to be taken 



3.2. The influence of Gould's Belt Star-Burst Region 
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Fig. 4. — Temporal variation of the 16 O CR primary density at the position of the Sun, for 
lOGeV per nucleon (upper panel) and 5TeV per nucleon (lower panel). 
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Fig. 5. — The density of 16 at E — 10 GeV per nucleoli in a 400 pc x 400 pc section of 
the Galactic plane (z = 0) during a local SN event for several times (in Myr). The Sun 
(r = 8.5 kpc, (j) = 1.025) is positioned at the center. The local CR density is strongly 
influenced by nearby SN, although the excess quickly disappears. 
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into account. We assume that the enhanced supernova rate applies to the entire volume of 
Gould's belt, i.e. we neglect a possible delay on account of the main-sequence phase of very 
massive stars. The probability distribution for supernova explosion would then be given by 
the time-dependent expansion model of Perrot & Grenier (2003) in addition to the stationary 
large-scale SN distribution in the Galaxy (see Eq.5). 

The results displayed in Fig. 6 show an enhanced variation of the local CR density due 
to the locally increased SN rate compared to the case without Gould's Belt, that has been 
given in the last section and is shown in Fig. 4. The locally enhanced SN rate also leads to 
a locally increased mean CR density which is most pronounced at low energies or, in other 
words, SNR need less efficiency as CR accelerators to replenish the galactic cosmic rays. 

3.3. Secondary Cosmic Ray Particles 

In section 3.1 we showed that the density of CR primary nuclei may vary up to about 
one order of magnitude in space and time. To investigate to what extent this is also true for 
secondary CR nuclei and, therefore, whether there are any effects on the ratio of secondary to 
primary CR isotopes, we further performed calculations including secondary nuclei. As the 
ratio of Boron to Carbon is widely used to quantify the parameters of various propagation 
models, we considered the isotopes 12 C, which was assumed to be pure primary and n B which 
was assumed to be produced solely by interactions of the primary 12 C with the interstellar 
matter. The results of these calculations are shown in Fig. 7, where we plot, as in Fig. 4, 
the CR density at 10 GeV per nucleon versus time for the primary isotope 12 C (upper panel) 
and the secondary isotope n B (lower panel), both at the position of the Sun. Fig. 7 shows 
that, although the density of the primary CR varies by a factor of 2, as expected from the 
previous calculations, the variation of the secondary CR particles is only a few per cent. 
This statement holds true also for higher particle energies. 

Having analyzed the variations of the CR density with time (Fig. 4, Fig. 7) and in the 
Galactic plane (Fig. 5), we now are interested in the variations of the CR flux perpendic- 
ular to the Galactic disk. In Fig. 8, we compare the densities of primary (upper panels) 
and secondary (lower panels) CR perpendicular to the Galactic plane, for CR energies of 
10 GeV (left panels) and 5 TeV (right panels), respectively, at four different instances of time, 
indicated by different line-styles. 

As in the case of Fig. 7, where we plotted the density of CR primary and secondary 
particles at one position in the Galactic plane versus time, these figures show, that there is 
only a marginal variation of the CR secondary distribution perpendicular to the Galactic 
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Fig. 6. — Same as Fig. 4, assuming a locally enhanced SN rate in Gould's Belt. 
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Fig. 7. — Temporal variation of the CR primary density 12 C (upper panel) and secondary 
n B (lower panel) at the position of the Sun at energy of 10 GeV per nucleon. The deviation 
from the average value for the CR secondary density is way smaller than that for the CR 
primary density. 
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Fig. 8. — The density of primary CRs (upper panel) and that of secondary CRs (lower panel) 
plotted as a function of distance from the Galactic plane at i?o = 8.5 kpc from the Galactic 
center. The left column is for a particle energy of 10 GeV per nucleon, and the right column 
applies to particles at 5 TeV energy. We show the density distributions at four arbitrarily 
chosen instances of time, indicated by different line-styles. 
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plane, despite rather large variations in the CR primary distributions. This behavior depends 
only weakly on the particle energy. Nevertheless, for particle energies of 5 TeV, there is a 
variation in the density of CR secondaries at z = of roughly 10%. These findings have 
far-reaching consequences for the ratio of secondary to primary CR, as discussed in the next 
section. 



4. Global signatures of an SNR origin of cosmic rays 

Equipped with a technique to calculate the spatial and temporal distribution of CR 
primary and secondary elements with SN-like objects as sources, we now discuss the obser- 
vational implications for CR astrophysics. We discuss the calculated spectra themselves and 
also compare them with measurements for the CR primary elements 16 O and 12 C. 

To obtain a measure of the variations of the CR spectra at the position of the Sun, 
we simulated 240 possible CR spectra both for the standard Galactic plane distribution of 
supernovae as given by Eq.5 and for the Gould's belt supernovae in addition to that. 

In Fig. 9 and Fig. 11 (upper panel), we give the range of possible spectra at the position of 
the Sun for the primary nuclei 16 and 12 C, respectively, for the case of randomly distributed 
sources (upper panels) and including also Gould's Belt (lower panels). In the same figures 
we compare the range of calculated spectra with data by Engelmann et al. (1990) (x), 
Muller et al. (1991) (□), Orth et al. (1978) (O) and Simon et al. (1980) (A). Those authors 
provide data for Carbon as well as oxygen. We model the effect of solar modulation using 
the force field approximation (Gleeson & Axford 1968), assuming a modulation parameter of 
$ = 500 MV. These locally measured data fit quite well in our calculated range of possible 
spectra. 

The amplitude of variations in the possible primary CR nuclei spectra is clearly seen and 
does not increase much with energy as is the case, e.g., for the electrons (section 4.3). The 
level of computational noise is very much less than the calculated fluctuation amplitudes 
for the primary CRs, as is indicated by the small level of fluctuations calculated for the 
secondary CRs in section 4.2. The local clustering of sources (section 4.1) due to Gould's 
Belt enhances the fluctuation amplitude, so the model fits the measurements even better. 
We also note a slight steepening of the averaged spectrum. 

Fig. 10 shows some examples of our calculated 12 C spectra. Note that the spectra mostly 
vary in their amplitudes, and only weakly in their forms. We now discuss these findings in 
detail. 
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Fig. 9. — The range of possible 16 spectra with (lower panel) and without (upper panel) 
Gould's Belt compared with measurements taken with balloons or satellites. The dark grey 
band shows the 68% containment probability range at each given energy and the light grey 
band gives the 95% range. The solid line marks the averaged (steady-state) spectrum. All 
spectra are as at the top of the atmosphere. 
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Fig. 11. — Same as Fig. 9 but for 12 C (upper panel) and n B, which is assumed to be purely 
secondary, produced from primary 12 C only (lower panel). 
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4.1. Effects due to the Clustering of Sources 

We investigated the effects of a local clustering of sources using Gould's Belt as an 
example. Our results are shown in Fig. 9 for the CR primary nucleus 16 and in Fig. 11 for 
CR secondary nucleus n B in each case with and without considering Gould's Belt. Local 
clustering of discrete CR sources leads to an enhanced fluctuation margin in the CR primary 
spectra and on average to a slightly steeper spectrum. The slight steepening of the spectrum 
is caused by the energy dependence of diffusion: the high-energy particles fill a somewhat 
larger volume than do their low-energy siblings. Then the differential number density must 
be softer than the total number spectrum, except for the short time interval the low-energy 
particle need to diffusively propagate to the location at which the spectrum is measured. 
Almost no effects of Gould's Belt can be seen in the CR secondary spectra. 

4.2. Implications on the Secondary to Primary Ratio 

Comparing the range of possible spectra for primary CR (see Fig. 9) and secondary 
CR (see Fig. 11) reveals that the variation in the possible CR nuclei spectra is much higher 
for primary nuclei than for secondary nuclei. Thus, the variation in the primary spectra 
should also be visible in the secondary-to-primary ratios. As the variation in the amplitude 
of the CR primary spectra is stronger than that of the spectral form, one should expect the 
secondary-to-primary ratios to vary mostly in amplitude but not so much in the spectral 
form. Taillet & Maurin (2003) show using a steady-state scenario that the CR measured 
at the earth only probe the CR propagation parameters in a rather small domain, and thus 
outside this domain the CR fluxes and the secondary-to-primary ratios may be different. 
Our calculations reveal that even if the propagation parameters do not vary, there is a 
variation in the secondary-to-primary ratio due to discrete sources. So provided CR nuclei 
are accelerated at SNR or other kinds of discrete sources, one must account for fluctuations 
in the ratio of secondary to primary CR particles, e.g. the Boron-to- Carbon ratio, which is 
widely used to determine the parameters of CR propagation models. In particular, the ratio 
should decrease in the vicinity of a source. Thus it is important to know whether or not we 
live in the vicinity of a recent supernova, as has been proposed by some authors (see Erlykin 
& Wolfendale (1997) and subsequent papers). 
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4.3. Comparison With Cosmic Ray Electrons 

We now compare our results derived for CR nuclei with the findings of Pohl & Esposito 
(1998) for CR electrons. Evidently, the fluctuation amplitude in the CR electron spectra 
strongly grows with energy, whereas that of CR nuclei spectra hardly depends on the particle 
energy (see e.g. Fig. 11). Kinks and dents are possible in CR electron spectra but are 
virtually not seen in the spectra of CR nuclei as shown in Fig. 10. 

5. Summary and discussion 

We have studied the impact of CR acceleration in SNR on the spectra of CR nuclei 
in the Galaxy. We found strong evidence that this assumption leads to CR spectra, which 
show significant variations in space and time. The behavior of the CR nuclei resembles 
that of protons, as suggested by first computations (Strong & Moskalenko 2001), but differs 
considerably from that of CR electrons (Pohl & Esposito 1998) . 

We have developed a method to numerically solve the propagation equation based on a 
series expansion, which allows us to use analytical solutions for part of the problem and an 
efficient distributed computing for the remainder. This method was employed to calculate 
the CR densities in the Galactic disk with high spatial (Az = 20 pc) and temporal (At = 
1 kyr) resolution for the primary nuclei 16 0, 12 C and the secondary nucleus n B for various 
distributions of SNR. This method allowed for the first time to calculate the CR densities 
in the Galactic disk with the high spatial and temporal resolution required to follow the 
CR density fluctuations caused by an SNR origin. We also studied the impact of a locally 
enhanced SN rate within Gould's Belt, a nearby star-forming region. 

We found strong variations of the CR nuclei flux by typically 20% with occasional 
spikes of much higher amplitude, but only minor changes in the spectral distribution. The 
locally measured primary CR spectra fit well into the obtained range of possible spectra. We 
further showed that the spectra of the secondary element Boron show almost no variations, 
so that the above findings also imply significant fluctuations of the Boron-to-Carbon ratio. 
Therefore the commonly used method of determining CR propagation parameters by fitting 
secondary-to-primary ratios appears flawed on account of the variations that these ratios 
would show throughout the Galaxy. 

Some indications that the CR flux varies in the Galaxy are given by the observation of 
the diffuse 7-ray emission in our Galaxy. Digel et al. (2001) performed 7-ray observations of 
the outer Galaxy. Their analysis of the diffuse 7-ray emission from giant molecular clouds in 
the Monoceros region suggests that the CR flux in the local Galactic arm and the neighboring 
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Perseus arm differ, i.e. the data suggests an enhancement of the CR density in the Perseus 
arm. These observations would be well explained by the variations in the flux of CR nuclei 
that results from a SNR origin. 

Hunter et al. (1997) find that the spectrum of the diffuse 7-ray emission toward the inner 
Galaxy can not be explained by the assumption that the locally observed CR spectrum and 
electron-to-proton ratio hold throughout the Galaxy. At 7-ray energies above some GeV, 
where the spectrum is presumably dominated by CR-nucleon-induced radiation, they find 
that the 7-ray flux measured exceeds by about 50% that expected if the local CR spectrum 
and electron-to-proton ratio hold throughout the Galaxy. These findings are difficult to 
explain by the calculations presented above as being caused by spectral variations of the CR 
nuclei flux, for we find these to be rather small. We have, however, neglected the possibility 
of a dispersion in the source spectral indices, which has been proposed as an explanation of 
the GeV excess (cf. Busching et al. 2001). 
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A. Green's Function for the Spherical Symmetric Propagation Equation 

In case of spherical symmetry, the propagation equation in spherical coordinates reads 
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(Eq. (15)): 




(Al) 



Using the ansatz 



N 




(A2) 



leads to an equation for = 0(r, ip, t) 




(A3) 
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For this equation, a Green's function can be found in the literature (Butkovskiy 1982). It 
reads 

g = e(t-t 



An^/rcy/k r r Q y/t - t 



Re-substituting Eq. (A2), we finally get 



xexp f- ,:' + r °M sinh f „ 2 r r ° (A5) 



4k (t- t ) J ' \4A: — t ) 

which we now show to be indeed the desired fundamental solution of Eq. (Al). For the 
derivatives we have: 

9_G = _1 r r ^(^)) c fA6) 

<9r r 2k(t-t ) 2k (t - 1 ) ginh f 2rr \ 1 ; 

<9r 2 r 2 + 2A;(t-t ) fc r (t - 1 ) ginh f 2rm \ ^\2k(t-t )J 

\4k(t-t ) J 

J i V cosh (jPg) „ , f r„ y 



^4fc(t-t ) 



cosh f 



1 c | r 2 + r 2 G 2rr cosn ^4fc(f-f ) / ^, 

^ 2 (* - *o) 4A; (t - t ) 2 4A; (f - t f sinh 

+ ^-^e(^)- (A8) 

Inserting the above into the Eq. (Al) yields 
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m sh ( 2rr ° 



2 ™sh[-^_l / r ^ 2 



2k (t- t ) J S inhf 2rr " ^ \2k(t-t ) 

4 G ._i_ fl+ _!L_ !^M G L W+S (A9) 
*(t-to) rMt-iolsinh^.^) J 

Comparing terms, we get 

r -i _zi±zjL sinh f 2rr o ^ , *(*-*o) | G _ 5 

\ 2 (t - t ) 4A: (t - t ) 2Sm \4k(t-t )J (t - t ) J 
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r 2 k(t-t ) \Ak(t-t ) 
So finally we have 



" e (( - t „) G - < A11 > 

Here S is the delta function for spherical polar coordinates, so we have to verify that the 
r.h.s. of Eq. (All) is a representation of the delta function 

S = ^5(r-r )5(t-t ). (A12) 
So with Eq. (A5), replacing the sinh by its definition, we have: 

S = 5(t-t ) 1 ^= exp (-&(*-*„)) (A13) 

(r-r ) 2 \ ( (r + r o y 

x exp — — exp 1 



4k (t- t ) J \ Ak(t-t ) 

Performing the limit for t — > t on the r.h.s. leads to 

lim - — X —= \== exp (-b (t - t )) 
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as required for a delta function. To check the normalization of Eq. (A14), we have to integrate 
over the spatial domain 



oo / 1 2-7r />n 



Gr 2 sin OdrdipdO 



r=0 Jtp=0 Je=o 
poo 

I = 4tt / Gr 2 dr 

Jr=0 



t=t 



t=t 
r 



xexp 



r 2 + rl 



sinh 



4k (t- t ) J \Ak(t-t ) 



2rr 



dr 



t=t 



With the integral 3.562.3 of Gradstein & Ryshik (1981) 

PCX) 

/ r exp (— f3 r) sinh (7 r) rfr = 

we finally arrive at 
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(A20) 
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which shows that the normalization is correct. 
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